% 旋转载体多天线中频信号仿真（基于文档1.2.2节准确相位公式）
% 核心：仅使用中频信号（无基带下变频），验证旋转相位与功率特性

%% 1. 加载配置（信号+接收机，确保相干积分时间1ms）
IF_carrierFreq = signalConfig.signal.CenterFrequency; % 中频载波频率(Hz)
sim_totalDuration = signalConfig.sim.duration;   % 总仿真时长
sim_startTime = 0;      % 开始时间
M = signalConfig.antennas.count;                  % 天线数量

%% 2. 初始化流处理模型（运动+信号特性，仅输出准确相位）
motionModel = MotionModel(signalConfig, rxConfig);       % 流处理运动模型（输出ψₘ(t)）
rotationModel = RotationSignalModel(signalConfig, rxConfig); % 输出准确相位


%% 3. 验证 motionModel 的输出结果(主要是卫星相对天线的俯仰角， 方位角的计算结果)并可视化，sim时间是 sim_totalDuration
resultCache = struct();
resultCache.timeTotal = [];                  % 总时间向量
resultCache.thetaMTotal = [];                  % 卫星相对天线的俯仰角（ψₘ(t)）
resultCache.phiMTotal = [];                  % 卫星相对天线的方位角（φₘ(t)）
resultCache.cosThetaMTotal = [];                  % 卫星相对天线的俯仰角（ψₘ(t)）的余弦值
resultCache.cosPhiMTotal = [];                  % 卫星相对天线的方位角（φₘ(t)）的余弦值
resultCache.phaseRotTotal = [];         % 旋转相位（rotationModel输出，文档1.2.2节公式）
resultCache.phaseApproxTotal = [];         % 近似相位（文档1.2.2节近似公式）
resultCache.powerTotal = [];            % 旋转对接收信号载波功率的影响（文档1.2.2节公式）
resultCache.EthetaTotal = [];            % 天线方向图 Etheta
resultCache.EphiTotal = [];            % 天线方向图 Ephi
for t =  sim_startTime : motionModel.config.motion.motionUpdateTime : (sim_startTime + sim_totalDuration)
    motionData = motionModel.step(t);
    rotationData = rotationModel.step(motionData);
    resultCache.timeTotal = [resultCache.timeTotal, t];
    resultCache.thetaMTotal = [resultCache.thetaMTotal, motionData.thetaM];
    resultCache.phiMTotal = [resultCache.phiMTotal, motionData.phiM];
    resultCache.cosThetaMTotal = [resultCache.cosThetaMTotal, cos(motionData.thetaM)];
    resultCache.cosPhiMTotal = [resultCache.cosPhiMTotal, cos(motionData.phiM)];
    resultCache.phaseRotTotal = [resultCache.phaseRotTotal, rotationData.phaseRot];
    % resultCache.phaseApproxTotal = [resultCache.phaseApproxTotal, rotationData.phaseApprox];
    resultCache.powerTotal = [resultCache.powerTotal, rotationData.powerRot];
    resultCache.EthetaTotal = [resultCache.EthetaTotal, rotationData.Etheta];
    resultCache.EphiTotal = [resultCache.EphiTotal, rotationData.Ephi];
end

% 卫星相对载体旋转中心的俯仰角（ψₘ(t)）和方位角（φₘ(t) 的变化）采用getElevationOrAzimuth函数
elevation = getElevationOrAzimuth(signalConfig.motion.elevMax, signalConfig.motion.elevMin, signalConfig.motion.passTime, resultCache.timeTotal, signalConfig.motion.useSineForElevation);  
azimuth = getElevationOrAzimuth(signalConfig.motion.azEnd, signalConfig.motion.azStart, signalConfig.motion.passTime, resultCache.timeTotal, signalConfig.motion.useSineForAzimuth);

figure;
subplot(2, 1, 1);
plot(resultCache.timeTotal, elevation);
title('卫星相对载体旋转中心的俯仰角（\theta_m(t)）');
xlabel('时间（秒）');
ylabel('俯仰角（弧度）');

% 卫星相对载体旋转中心的方位角（φₘ(t) 的变化）
subplot(2, 1, 2);
plot(resultCache.timeTotal, azimuth);
title('卫星相对载体旋转中心的方位角（\phi_m(t)）');
xlabel('时间（秒）');
ylabel('方位角（弧度）');


% 卫星相对天线相位中心的俯仰角（ψₘ(t)）和方位角（φₘ(t) 的变化）
figure;
subplot(2, 1, 1);
plot(resultCache.timeTotal, resultCache.thetaMTotal);
title('卫星相对天线相位中心的俯仰角（\theta_m(t)）');
xlabel('时间（秒）');
ylabel('俯仰角（弧度）');
legend(arrayfun(@(x) sprintf('天线 %d', x), 1:M, 'UniformOutput', false));

subplot(2, 1, 2);
plot(resultCache.timeTotal, resultCache.phiMTotal);
title('卫星相对天线相位中心的方位角（\phi_m(t)）');
xlabel('时间（秒）');   
ylabel('方位角（弧度）');
legend(arrayfun(@(x) sprintf('天线 %d', x), 1:M, 'UniformOutput', false));

% 卫星相对天线相位中心的俯仰角（ψₘ(t)）和方位角（φₘ(t) 的余弦值）       
figure;
subplot(2, 1, 1);
plot(resultCache.timeTotal, resultCache.cosThetaMTotal);
title('卫星相对天线相位中心的俯仰角（\theta_m(t)）的余弦值');
xlabel('时间（秒）');
ylabel('俯仰角（弧度）的余弦值');
legend
legend(arrayfun(@(x) sprintf('天线 %d', x), 1:M, 'UniformOutput', false));

subplot(2, 1, 2);
plot(resultCache.timeTotal, resultCache.cosPhiMTotal);
title('卫星相对天线相位中心的方位角（\phi_m(t)）的余弦值');
xlabel('时间（秒）');   
ylabel('方位角（弧度）的余弦值');
legend(arrayfun(@(x) sprintf('天线 %d', x), 1:M, 'UniformOutput', false));          

% 卫星相对天线相位中心的俯仰角（ψₘ(t)）和方位角（φₘ(t) 的正弦值）
figure;
subplot(2, 1, 1);
plot(resultCache.timeTotal, sin(resultCache.thetaMTotal));
title('卫星相对天线相位中心的俯仰角（\theta_m(t)）的正弦值');
xlabel('时间（秒）');
ylabel('俯仰角（弧度）的正弦值');
legend(arrayfun(@(x) sprintf('天线 %d', x), 1:M, 'UniformOutput', false));

subplot(2, 1, 2);
plot(resultCache.timeTotal, sin(resultCache.phiMTotal));
title('卫星相对天线相位中心的方位角（\phi_m(t)）的正弦值');
xlabel('时间（秒）');   
ylabel('方位角（弧度）的正弦值');
legend(arrayfun(@(x) sprintf('天线 %d', x), 1:M, 'UniformOutput', false));          


% 可视化：旋转对接收信号载波功率的影响（P(t)）
figure;
subplot(2, 2, 1);
plot(resultCache.timeTotal, resultCache.phaseRotTotal);
title('旋转对接收信号载波相位的影响（\phi_{rot}(t)）');
xlabel('时间（秒）');
ylabel('相位（弧度）');
legend(arrayfun(@(x) sprintf('天线 %d', x), 1:M, 'UniformOutput', false));
subplot(2, 2, 2);
plot(resultCache.timeTotal, resultCache.powerTotal);
title('旋转对接收信号载波功率的影响（P(t)）');
xlabel('时间（秒）');
ylabel('贴片天线方向图');
legend(arrayfun(@(x) sprintf('天线 %d', x), 1:M, 'UniformOutput', false));

% 补充可视化: 求旋转对接收信号载波相位的影响（\phi_{rot}(t)）的幅度谱 以及 载波幅度的幅度谱
subplot(2, 2, 3);
plotSingleSidedSpectrum(resultCache.phaseRotTotal, 1 / motionModel.config.motion.motionUpdateTime);
title('\phi_{rot}(t)的幅度谱');
subplot(2, 2, 4);
plotSingleSidedSpectrum(sqrt(resultCache.powerTotal), 1 / motionModel.config.motion.motionUpdateTime);
title('P(t)的幅度谱');


%% 功率 + 俯仰角双纵轴图
figure('Color','w', 'Position',[200 200 860 420]);

% 左轴：归一化接收功率（建议画成 dB，更直观）
yyaxis left                                      
plot(resultCache.timeTotal, ...
     10*log10(resultCache.powerTotal), ...   % 直接转 dB
     'LineWidth', 1.6);
ylabel('接收功率 (dB)','FontWeight','bold');
ylim([10*log10(min(resultCache.powerTotal(:)))-2 , 2]);  % 自动留点余量
grid on;

% 右轴：卫星俯仰角（所有天线基本一致，只画第1根就够，多了乱）
yyaxis right
elev_deg = rad2deg( resultCache.thetaMTotal(1,:) );   % 取第1根天线俯仰角即可
plot(resultCache.timeTotal, elev_deg, 'k--', 'LineWidth', 2);
ylabel('卫星俯仰角 Elevation (deg)','FontWeight','bold','Color','k');
ylim([0 95]);

% 标题、坐标轴、图例
title('旋转过程中天线接收功率与卫星俯仰角变化关系','FontSize',13,'FontWeight','bold');
xlabel('时间 (s)');
legend('归一化接收功率 P(t)', '卫星俯仰角 el(t)', ...
       'Location','southoutside','FontSize',11);

       
% 补充可视化，天线图 Etheta 和 Ephi 以及 天线方向图 Etheta 和 Ephi 的幅度谱
figure;
subplot(2, 2, 1);
plot(resultCache.timeTotal, resultCache.EthetaTotal);
title('天线方向图 Etheta');
xlabel('时间（秒）');
ylabel('Etheta');
legend(arrayfun(@(x) sprintf('天线 %d', x), 1:M, 'UniformOutput', false));
subplot(2, 2, 2);
plotSingleSidedSpectrum(resultCache.EthetaTotal, 1 / motionModel.config.motion.motionUpdateTime);
title('Etheta 的幅度谱');

subplot(2, 2, 3);
plot(resultCache.timeTotal, resultCache.EphiTotal);
title('天线方向图 Ephi');
xlabel('时间（秒）');
ylabel('Ephi');
legend(arrayfun(@(x) sprintf('天线 %d', x), 1:M, 'UniformOutput', false));
subplot(2, 2, 4);
plotSingleSidedSpectrum(resultCache.EphiTotal, 1 / motionModel.config.motion.motionUpdateTime);
title('Ephi 的幅度谱');

%% 补充可视化: 方向图
ref_ant = 1;

% 第1步：取出俯仰角（rad）
elev_rad = resultCache.thetaMTotal(ref_ant, :);        % 俯仰角，90°=天顶

% 第2步：转换为天顶角
zenith_rad = elev_rad;                          % 天顶角，0°=天顶

% 第3步：对应的实际接收功率（已经与时间严格对齐）
power_lin = resultCache.powerTotal(ref_ant, :);
power_dB  = 10*log10(power_lin);
power_dB  = power_dB - max(power_dB);                  % 归一化到 0 dB

% 第4步：理论方向图（经典圆极化贴片天线模型）
theta_model = 0 : 0.001 : pi/2;                        % 天顶角 0~90°
pattern_dB  = 10*log10( cos(theta_model).^3 );         % 功率方向图 G(θ) = cos³(θ)
pattern_dB  = pattern_dB - max(pattern_dB);

% 第5步：绘图
figure('Color','w', 'Position',[100 100 820 720]);
ax = polaraxes; hold on; grid on;

polarplot(theta_model, pattern_dB,  'k-',  'LineWidth', 2);     % 理论曲线
polarplot(zenith_rad,  power_dB,   'r-',  'LineWidth', 3);     % 实际轨迹（完美贴合！）

% 专业极坐标设置（缺一不可）
ax.ThetaZeroLocation = 'top';           % 0° 在正上方 = 天顶
ax.ThetaDir          = 'clockwise';     % 顺时针 = 天顶角增大
ax.ThetaLim          = [0 90];
ax.ThetaTick         = 0:15:90;
ax.ThetaTickLabel    = {'0°(天顶)','15°','30°','45°','60°','75°','90°(地平)'};
ax.RLim              = [-30 0];
ax.RTick             = -30:10:0;

title(sprintf('天线 %d：卫星实时天顶角在归一化方向图上的投影', ref_ant), ...
      'FontSize',14,'FontWeight','bold');

legend('理论方向图 G(θ_{zenith}) = cos^3(θ_{zenith})', ...
       '实际接收功率轨迹', ...
       'Location','southoutside', 'FontSize',12);

% 可选：标出 -3dB 点
% theta_3dB = acos( (0.5)^(1/3) );  % ≈ 65.4°
% polarplot([0 theta_3dB], [-3 -3], '--', 'Color',[0.5 0.5 0.5], 'LineWidth',1.2);
% text(theta_3dB+3, -3.5, ' -3 dB', 'FontSize',11, 'FontWeight','bold');